      subroutine parrefl_diffuse(ncellsp, ijkctmp, iphead, npart, itdim, nsampl&
         , vrms, t_wall, ico, iwid, jwid, kwid, x, y, z, wate, bxv, byv, bzv, &
         px, py, pz, up, vp, wp, pxi, peta, pzta) 
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use modify_com_M 
 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: ncellsp 
      integer , intent(in) :: npart 
      integer  :: itdim 
      integer , intent(in) :: nsampl 
      integer , intent(in) :: iwid 
      integer , intent(in) :: jwid 
      integer , intent(in) :: kwid 
      integer , intent(in) :: ijkctmp(*) 
      integer , intent(in) :: iphead(*) 
      integer , intent(in) :: ico(0:npart) 
      real(double) , intent(in) :: vrms(*) 
      real(double) , intent(in) :: t_wall(*) 
      real(double) , intent(in) :: x(*) 
      real(double) , intent(in) :: y(*) 
      real(double) , intent(in) :: z(*) 
      real(double)  :: wate(itdim,*) 
      real(double)  :: bxv(*) 
      real(double)  :: byv(*) 
      real(double)  :: bzv(*) 
      real(double) , intent(inout) :: px(0:npart) 
      real(double) , intent(inout) :: py(0:npart) 
      real(double) , intent(inout) :: pz(0:npart) 
      real(double) , intent(inout) :: up(0:npart) 
      real(double) , intent(inout) :: vp(0:npart) 
      real(double) , intent(inout) :: wp(0:npart) 
      real(double) , intent(in) :: pxi(0:npart) 
      real(double) , intent(in) :: peta(0:npart) 
      real(double) , intent(in) :: pzta(0:npart) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: n, ijk_host, np, inew, jnew, knew, ijk_chek, nbin 
      real(double) :: normx, normy, normz, pi, tsi, eta, dx1, dx2, dy1, dy2, &
         dz1, dz2, x_intersect, y_intersect, z_intersect, eta_r, tsi_r, &
         x_random, y_random, z_random, rnorm, dx, dy, dz, dxdotn, udotn_old, th&
         , fi, vmag, udotn 
!-----------------------------------------------
!   E x t e r n a l   F u n c t i o n s
!-----------------------------------------------
      real(double) , external :: ranf 
!-----------------------------------------------
!
!     ***************************************************************
!
!     a routine for diffuse reflection of particles
!     from a wall at temperature T_wall
!     particle is created with coordinates of intersection
!     with wall
!     new velocity calculated from Maxwellian distribution
!
!     **************************************************************
!
!
!
      pi = acos(-1.) 
!
!dir$ ivdep
!
      do n = 1, ncellsp 
!
         ijk_host = ijkctmp(n) 
         np = iphead(ijk_host) 
!
         inew = int(pxi(np)) 
         jnew = int(peta(np)) 
         knew = int(pzta(np)) 
!
         ijk_chek = (knew - 1)*kwid + (jnew - 1)*jwid + (inew - 1)*iwid + 1 
!
         tsi = pxi(np) - inew 
         eta = peta(np) - jnew 
!
         dx1 = (1. - eta)*(x(ijk_chek+iwid+kwid)-x(ijk_chek+kwid)) + eta*(x(&
            ijk_chek+iwid+jwid+kwid)-x(ijk_chek+jwid+kwid)) 
!
         dx2 = (1. - tsi)*(x(ijk_chek+jwid+kwid)-x(ijk_chek+kwid)) + tsi*(x(&
            ijk_chek+iwid+jwid+kwid)-x(ijk_chek+iwid+kwid)) 
!
         dy1 = (1. - eta)*(y(ijk_chek+iwid+kwid)-y(ijk_chek+kwid)) + eta*(y(&
            ijk_chek+iwid+jwid+kwid)-y(ijk_chek+jwid+kwid)) 
!
         dy2 = (1. - tsi)*(y(ijk_chek+jwid+kwid)-y(ijk_chek+kwid)) + tsi*(y(&
            ijk_chek+iwid+jwid+kwid)-y(ijk_chek+iwid+kwid)) 
!
         dz1 = (1. - eta)*(z(ijk_chek+iwid+kwid)-z(ijk_chek+kwid)) + eta*(z(&
            ijk_chek+iwid+jwid+kwid)-z(ijk_chek+jwid+kwid)) 
!
         dz2 = (1. - tsi)*(z(ijk_chek+jwid+kwid)-z(ijk_chek+kwid)) + tsi*(z(&
            ijk_chek+iwid+jwid+kwid)-z(ijk_chek+iwid+kwid)) 
!
         normx = dy1*dz2 - dy2*dz1 
         normy = dz1*dx2 - dz2*dx1 
         normz = dx1*dy2 - dx2*dy1 
!
         x_intersect = tsi*(x(ijk_chek+iwid+kwid)*(1.-eta)+x(ijk_chek+iwid+jwid&
            +kwid)*eta) + (1. - tsi)*(x(ijk_chek+kwid)*(1.-eta)+x(ijk_chek+jwid&
            +kwid)*eta) 
!
         y_intersect = tsi*(y(ijk_chek+iwid+kwid)*(1.-eta)+y(ijk_chek+iwid+jwid&
            +kwid)*eta) + (1. - tsi)*(y(ijk_chek+kwid)*(1.-eta)+y(ijk_chek+jwid&
            +kwid)*eta) 
!
         z_intersect = tsi*(z(ijk_chek+iwid+kwid)*(1.-eta)+z(ijk_chek+iwid+jwid&
            +kwid)*eta) + (1. - tsi)*(z(ijk_chek+kwid)*(1.-eta)+z(ijk_chek+jwid&
            +kwid)*eta) 
!
         eta_r = ranf() 
         tsi_r = ranf() 
!
         x_random = tsi_r*(x(ijk_chek+iwid+kwid)*(1.-eta_r)+x(ijk_chek+iwid+&
            jwid+kwid)*eta_r) + (1. - tsi_r)*(x(ijk_chek+kwid)*(1.-eta_r)+x(&
            ijk_chek+jwid+kwid)*eta_r) 
!
         y_random = tsi_r*(y(ijk_chek+iwid+kwid)*(1.-eta_r)+y(ijk_chek+iwid+&
            jwid+kwid)*eta_r) + (1. - tsi_r)*(y(ijk_chek+kwid)*(1.-eta_r)+y(&
            ijk_chek+jwid+kwid)*eta_r) 
!
         z_random = tsi_r*(z(ijk_chek+iwid+kwid)*(1.-eta_r)+z(ijk_chek+iwid+&
            jwid+kwid)*eta_r) + (1. - tsi_r)*(z(ijk_chek+kwid)*(1.-eta_r)+z(&
            ijk_chek+jwid+kwid)*eta_r) 
!
!
         rnorm = 1./(normx**2 + normy**2 + normz**2) 
!
         dx = px(np) - x_intersect 
         dy = py(np) - y_intersect 
         dz = pz(np) - z_intersect 
!
         dxdotn = (dx*normx + dy*normy + dz*normz)*rnorm 
!
!     reflect particle inside
!
         px(np) = px(np) - 2.*normx*dxdotn 
         py(np) = py(np) - 2.*normy*dxdotn 
         pz(np) = pz(np) - 2.*normz*dxdotn 
!
!     shift to random point of origin on surface
!
         px(np) = px(np) - x_intersect + x_random 
         py(np) = py(np) - y_intersect + y_random 
         pz(np) = pz(np) - z_intersect + z_random 
         udotn_old = (up(np)*normx+vp(np)*normy+wp(np)*normz)*rnorm 
!     generate new particle velocities from Maxwellian
!     distribution
!
         th = 2.*pi*ranf() 
         fi = acos(2.*ranf() - 1.) 
         nbin = ifix(ranf()*float(nsampl - 1) + 1.) 
         vmag = vrms(nbin) 
!
         up(np) = t_wall(ico(np))*vmag*cos(th)*sin(fi) 
         vp(np) = t_wall(ico(np))*vmag*sin(th)*sin(fi) 
         wp(np) = t_wall(ico(np))*vmag*cos(fi) 
!
         udotn = (up(np)*normx+vp(np)*normy+wp(np)*normz)*rnorm 
!
         if (udotn <= 0.) cycle  
!
         up(np) = up(np) - 2.*udotn*normx 
         vp(np) = vp(np) - 2.*udotn*normy 
         wp(np) = wp(np) - 2.*udotn*normz 
!
      end do 
!
      return  
      end subroutine parrefl_diffuse 
